function sa2(N)

  

  
  time_steps = [0:N-1]./(N-1)*0.999;

  function resSa2b = Sa2b(z)
    if (z < 0.0001)
      resSa2b = (-1/3)+(3/8).*z+(-5/12).*z.^3+(3/7).*z.^4+(-4/9).*z.^6+(9/20).*z.^7+(-11/24).*z.^9+(6/13).*z.^10+(-7/15).*z.^12+(15/32).*z.^13+(-17/36).*z.^15+(9/19).*z.^16+(-10/21).*z.^18+(21/44).*z.^19;
    else
      resSa2b = z.^(-2).*((1/2)+(-1/12).*z.^(-1).*(3.^(1/2).*pi+6.*z.^3.*(1+z+z.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)));
    end


  end
  function resM2b = M2b(z)
  % pulled out z^3
    resM2b = zeros(2,2);
    resM2b(1,1) = 1;
    %resM(2,1) = 2 + z.^3;
    %resM(2,2) = z.*(z.^3 - 1);
    resM2b(2,1) = (-4+7.0.*z^3);
    resM2b(2,2) = z.*(z.^3-1);
  end
  function res2b = deq2b(z, u)
    a   = u(1);
    az  = u(2);
    res2b = zeros(2,1);
    res2b(1) = az;
    res2b(2) = 2*Sa2b(z) - 9.0.*z.^2.*a + 2.0.*z;
  end
  %uin = [(2+pi./3./sqrt(3)-log(3))./3; -4./6];
  uin = [(1/3)+(1/108).*((-1).*3.^(1/2).*pi+9.*((-2)+log(3))); -1./4];
  uin
  opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @M2b);
  [tout, uout] = ode15s(@deq2b, time_steps, uin, opts);
  %plot(tout, uout(:,1)./log(1-tout))
  plot(uout)

  %uout(:,1)./log(1-tout)


end
